【GIS实验课】05 坡度提取的不确定性分析（分享批量重采样、批量掩膜、坡度批量提取代码）

今日分享：

后台回复“批量”可以获取批量重采样、批量掩膜、批量坡度提取和批量分区统计的代码，不过你们懂得。

01主要内容

本次实验下载的是GDEMV2 30M分辨率数字高程数据，利用Python提取不同分辨率的DEM，基于上述不同分辨率DEM提取每种地貌类型的平均坡度，最后以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图，并对散点图进行拟合，通过回归算法求得回归方程的系数及常数项。

02具体要求

1.以30m空间分辨率的DEM数据为基础数据，重采样为40、50、60、70、80、90、100、110、120 m共10组不同分辨率的DEM。

2. 基于不同分辨率DEM提取每种地貌类型的平均坡度。

3. 以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图，并对散点图进行拟合，通过回归算法求得回归方程的系数及常数项。

03具体步骤
◐ 1. 使用ArcPy进行处理

1.1 将五景DEM数据镶嵌起来然后利用ArcPy进行批量重采样，具体代码如下所示：

import arcpyin_raster = r"C:\Users\Admin\Desktop\GIS Practice\ dem_mosaic.tif"out_raster_workspace = "C:\\Users\\Admin\\Desktop\\GISPractice\\ resample\\"for n in range(40,130,10): out_raster =out_raster_workspace + "resample_" + str(n) + ".tif" arcpy.Resample_management(in_raster, out_raster,str(n),"CUBIC")print "OK!"

1.2 将重采样得到10组不同分辨率的DEM，利用行政区的矢量边界，编写Python代码进行批量剪裁，具体代码如下所示：

import arcpy,os,globfrom arcpy import envarcpy.CheckOutExtension("Spatial")filepath=r"C:\\Users\\Admin\\Desktop\\GISPractice\\ resample"env.workspace = filepathos.chdir(filepath)Rasters = glob.glob("*.tif")outfile="C:\\Users\\Admin\\Desktop\\GISPractice\\test2\\clip"inMaskData = "C:\\Users\\Admin\\Desktop\\GISPractice\\hefei\\He.shp"for Raster in Rasters: inRaster =Raster arcpy.gp.ExtractByMask_sa(Raster, inMaskData, outfile+"\\"+Raster.split(".")[0]+'_clip.tif')print "ok"

批量剪裁结果：

 图1｜批量剪裁结果

1.3 将上述批量剪裁完的不同分辨率的DEM数据进行批量提取坡度，具体的Python代码如下所示:

import arcpyfrom arcpy import envenv.workspace = "C:\\Users\\Admin\\Desktop\\GISPractice\\ clip" Rasters =arcpy.ListRasters("*","tif")arcpy.CheckOutExtension("3d")for inRaster in Rasters: outMeasurement = "DEGREE" zFactor =1 outRaster = "C:\\Users\\Admin\\Desktop\\GISPractice\\ slope\\" +\ inRaster.split(".")[0] +"_Slope.tif" arcpy.Slope_3d(inRaster,outRaster, "DEGREE", zFactor)print "you are very wise!"

1.4 基于上述不同分辨率DEM提取每种地貌类型的平均坡度（批量分区统计），具体代码如下所示：

import arcpyfrom arcpy import envfrom arcpy.sa import *env.workspace ="C:\\Users\\Admin\\Desktop\\GIS Practice\\ slope"ValueRaster=arcpy.ListRasters("*","tif")arcpy.CheckOutExtension("Spatial")for inValueRasterin ValueRaster: inZoneData="C:\\Users\Admin\\Desktop\\GISPractice\\quanguodimao__WGS1984\\hefeidimao_prj.shp" zoneField = "GEOMOR_ID" outTable="C:\\Users\Admin\\Desktop\\GISPractice\\ZonalSta\\"+inValueRaster.split(".")[0]+".dbf" outZSaT=ZonalStatisticsAsTable(inZoneData,zoneField, inValueRaster,outTable,"NODATA","MEAN")print "you arevery wise!"

1.5 利用EXCEL统计不同分辨率DEM下提取的每种地貌类型的平均坡度，如下表所示：

表1｜不同分辨率DEM下提取的每种地貌类型的平均坡度

以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图，并对散点图进行拟合，通过回归算法求得回归方程的系数及常数项（使用的工具是excel）。

图2｜不同地貌类型的散点图

表2｜不同地貌拟合函数表达式

由图2和表2可以看出平均坡度值与分辨率呈现很强的线性关系且是正相关关系，可以以线性方程来描述这种关系方程式为y = ax + b，其中x为DEM分辨率，y为不同DEM分辨率下的平均坡度。低海拔丘陵的斜率最大，低海拔洪积平原的斜率最小，斜率绝对值之差为0.3888。从整体上看，按照拟合曲线的斜率，可大致将上述地貌类型分为两类：(1)斜率较大类：低海拔丘陵、低海拔冲积洪积台地、低海拔冲积平原、低海拔冲积扇平原；(2)斜率略小类：低海拔小起伏山地、低海拔冲积台地、低海拔洪积平原、低海拔湖积平原、湖泊。对某市各地貌类型拟合曲线的斜率求均值作为该研究区的平均曲率，求得a = -0.32493。

◐ 2. 用Model Builder进行处理

具体模型如下所示：

 图3｜模型示意

在Model Builder中拖入各种数据进行建模，先加入包含不同分辨率DEM数据的文件夹clip，然后插入栅格迭代器，并设置工作空间或栅格目录为带有迭代号的文件夹clip，接着加入按掩模提取工具，将某市区域提取出来，然后加入Slope工具和分区统计工具，在分区统计工具设置中，输入要素区域数据为某市地貌矢量数据，使用地貌数据的ID字段对每种分辨率下的坡度数据进行统计，输出文件的名称为：%名称%.tif，统计类型为均值，运行该模型,则可以得到相应数据类型下的平均坡度。

Tips：

在编写ArcPy代码进行DEM数据的批量重采样的时候出现了报错，经过排查发现主要原因是因为out_raster = out_raster_workspace +"resample_" + str(n) + ".tif"这一句代码出现了错误，我们对DEM数据进行重采样，从30米到120米一共有10景DEM数据，输出的每个DEM的名称肯定是不一样的，都是根据DEM数据的分辨率来进行命名，采用的Python语句是：for n in range(40,130,10)，而问题就是出现这里，这里面的n是表示数字，所以在下面的代码中需要写成str(n)，因为如果不这样写的话，这个n会被认定为一个无效字符。

除此之外，在利用矢量边界对不同分辨率的DEM进行批量剪裁的时候出现了错误，在这之前我也编写ArcPy做过不少批量剪裁，不过是用不同的矢量边界去裁剪同一个栅格，遍历矢量数据的语法是：Features=arcpy.ListFiles(“*shp”)，但是本次需要的是用同一个矢量边界去批量剪裁多个栅格数据，所以遍历数据的语句则改为: Rasters =glob.glob("*.tif")，在编写代码的时候我导入的库有：arcpy、os、glob，如果只是导入arcpy则程序无法执行，通过多次的调试代码终于运行成功!!!

后台回复“批量”可以获取批量重采样、批量掩膜、坡度批量提取和批量分区统计的代码，emmmmmm,不过你们懂得==

作者｜不许人间见白头

排版｜Moon

校阅｜数读菌、不许人间见白头

那今天就到这里结束啦，欢迎留言讨论。文中的图片文字未经许可不要随便“引用”。

如果可以的话，希望能够转发分享，点个在看并且点个赞，给个赞赏~~也欢迎规范转载~

也希望大家和我多留言互动啊！（据说这样可以增加我的推送在你的订阅号里出现的概率）

文章推荐 【数据分享】中国生态服务功能重要性与中国生态系统敏感性【双评价学习笔记】农业：环境评价【GIS实验课】04 GIS数据的批处理方法【GIS实验课】03 城镇住房选址适宜区评价分析【ENVI实验课】02利用ENVI进行农业耕作用地变化监测【ENVI实验课】01基于像元二分模型的植被覆盖度反演【GeoDa基础】GeoDa平台下的江西省县级统计数据分析【GIS实验课】02 超市选址适宜性评价【数据整理】2019年行政区划调整的数据更新与分享【双评价学习笔记】农业：气候评价（协同克里金方法）【GIS实验课】01山区土壤侵蚀评价及其地形的关系——以皖西大别山为例【GIS基础】基本空间分析工具：缓冲区、网络分析、相交分析、密度分析

▼

需要你的“分享”和“在看”

END>

如需全文转载文章、投稿或者合作

可添加微信

(回复超慢！！！)

（不要添加我问各种问题，我大概率不会的==）

（入群请一定要备注入群）

（添加后会在晚上非工作时间通过，请稍安勿躁）

公众号

微博

▼ 点击阅读原文，使用关键词搜索历史文章